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SIGNIFICANT ACCOMPLISHMENTS 


This Semiannual Report covers research conducted under NASA Contract NAG5-459 for the 
period 1 August 1984 through December 31, 1985. 

Our effort on behalf of the Crustal Dynamics Project during that period focused on the 
development of methodologies suitable for the analysis of space-geodetic data sets for the 
estimation of crustal motions, in conjunction with results derived from land-based geodetic data, 
neo-tectonic studies, and other geophysical data. These methodologies were used to provide 
estimates of both global plate motions and intraplate deformation in the western U.S. Significant 
accomplishments include: 

1. Results from the satellite ranging experiment for the rate of change of the baseline length 
between San Diego and Quincy, California indicated that relative motion between the North 
American and Pacific plates over the course of the observing period during 1972-1982 were 
consistent with esimates calculated from geologic data averaged over the past few million years. 

2. This result, when combined with other kinematic constraints on western U.S. 
deformation derived from land-based geodesy, neo-tectonic studies, and other geophysical data, 
places limits on the possible extension of the Basin and Range province, and implies significant 
deformation is occurring west of the San Andreas fault 

3. A new methodology was developed to analyze vector-position space-geodetic data to 
provide estimates of relative vector motions of the observing sites. The algorithm is suitable for the 
reduction of large, inhomogeneous data sets, and takes into account the full position covariances, 
errors due to poorly resolved earth orientation parameters and vertical positions, and reduces baises 
due to inhomogeneous sampling of the data. 

4. This methodology was applied to the problem of estimating the rate-scaling parameter of a 
global plate tectonic model using satellite laser ranging observations over a five-year interval. The 
results indicate that the mean rate of global plate motions for that interval are consistent with those 
averaged over several million years, and are not consistent with quiescent or greatly accelerated 
plate motions. 

5. This methodology was also used to provide constraints on deformation in the western 
U.S. using very long baseline interferometry observations over a two-year period. Motions 
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relative to a fixed North American plate frame were consistent with global plate tectonic models at 
Vandenburg, but were less at Monument Peak, Calif onia, suggesting significant deformation west 
of the San Andreas may be occurring between those sites. 


PROBLEMS AND RECOMMENDATIONS 
None 

DATA UTILITY 

Not applicable 

FUNDS EXPENDED 

As of 31 December 1985, a total of $72,720 had been spei^ out of the current fund 
limitation of $ 1 36,500. 
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APPENDIX I 


relative movement of tectonic plates 

IN CALIFORNIA OBSERVED B V 
SATELLITE LASER RANGING 


Abstract: A satellite laser ranging experiment conducted by NASA since 1972 
has measured the relative motion between the North America and Pacific plates in 
California. Based on these measurements, the 896-km distance between San Diego 
and Quincy, California, is shortening at 62 + 9 mm/yr. This geodetic estimate is 
consistent with the rate of motion between the two plates, calculated from geolog- 
ical data to be 53 + 3 mm/yr, averaged over the past few million years. 


Lasers capable of tracking near-earth orbiting satellites have been providing 
important geodetic data for mare than a decade. Beginning in 1972 with the San 
Andreas Fault Experiment (SAFE}, short-term tectonic motions along a baseline 
spanning the San Andreas system have been monitored by satellite laser ranging 
(SLR) in an effort to improve our understanding of earthquake hazards in California 
(Fig. 1). This experiment, which continues as part of NASA's Crustal Dynamics 
Project, has measured the variations in the 896-km distance between two tracking 
systems located on opposite sides of the fault with an accuracy unachievable by 
classical geodetic techniques. 

Previous reports (1, 2. 3) have shown the technique to be capable of monitor- 
ing short-term C*1Q yr) tectonic movements at the centimeter-per-year level. In 
this paper and the one which follows, we analyze the entire data set collected 
over the 11-year period 1972-1982 and discuss some of its tectonic implications. In 
particular, we shall show that space geodetic techniques are beginning to place 
more stringent constraints on large-scale deformation of the western United States 
than those derived from classical geological and geophysical observations (4). 

Mobile laser systems were used throughout the tracking campaigns, campaigns. 
Both sites were reoccupied approximately every two years for a duration of several 
months in order to satisfy certain criteria of data acquisition and, earlier on, 
tracking geometry. The first results were obtained by tracking a low orbiting 
spacecraft. Beacon Explorer-C (BE C), launched in 1965 and equipped with an array 
of laser retro-reflectors to provide, among other objectives, a spaceborne target for 
laser engineering studies. In 1976, a dedicated Laser Geodynamics Satellite 
(LAGEOS) was launched into a high earth orbit with the sole purpose of providing 
an optimal target for global laser tracking in the context of geodynamic investiga- ’ ^ 
tions. 


- 5 - 



Relative Movement of Tectonic Plates from SLR 


In satellite laser ranging, the observed quantities are the nearly instantaneous 
ranges from a ground station to a spacecraft and the times at which these dis- 
tances are measured. If the evolution of the orbit and the variations of the earth 
rotational parameters are known accurately, then the positions of the observing 
stations (to wi thin an arbitrary shift in absolute longitude) can be calculated in the 
geocentric reference frame used to describe the satellite motion. Analysis of 
satellite laser ranging data requires careful consideration of: \ 

(1) the accuracy with which the motion of a near-earth satellite subjected to 
complex gravitational and non-conservative forces can be modeled, 

(2) the accuracy of the observations themselves, which limits the available 
approaches, and 

(3) the requirements of globally distributed tracking support. 

BE-C, a low-altitude spacecraft with a complex shape, posed a difficult tra- 
jectory modeling problem. Orbital accuracy was severely degraded if more than a 
few revolutions were analyzed simultaneously; that is, analysis had to be restricted 
to relatively short arcs. A special data reduction technique was consequently dev- 
ised to accommodate data collected from two tracking sites only. In spite of its 
intrinsic limitations for measuring absolute baselines, this technique is capable of 
detecting changes in the inter-station distance by using uniform tracking geometry 
and identical force models for each bi-yeariy analysis (2, 5). 

The launch of LAGEOS at high altitude provided a well-designed laser target 
far less sensitive to pooriy known perturbing forces such as short wavelength 
gravity and atmospheric drag. With the deployment of a worldwide network of 
more advanced laser systems, trajectory accuracies achievable for LAGEOS 
approach the quality of the laser tracking itself, leading to improved estimates of 
the gravity field at long wavelengths and to highly accurate earth rotational 

parameters (6, 7). Furthermore, the SAFE line has been incorporated as part of 
the several hundred lines measured annually in the context of NASA's Crustal 

Dynamics Project, so that a special data reduction technique is no longer required. 
Unlike the baselines obtained from BE-C tracking, those derived from LAGEOS are 
believed to be accurate in an absolute sense (they do, however, depend on the 

values of GM^ and the speed of light used in the reduction procedure). 

LAGEOS observations are analyzed in various ways. We have derived robust 
three-dimensional station positions by combining twelve monthly orbital arcs into 
annual solutions. In this case, errors arising from weaker data sets and force 

model imperfections are reduced through averaging over an entire year. 

In addition to these annual solutions, we also computed selected subset solu- 
tions with finer temporal sampling, which extend the capabilities of SLR to resolve 
relatively short (<1000 km) baselines, using vastly smaller data sets ( from a few 
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weeks or less). These rely on a short-arc technique, called Baseline Estimation 
from Simultaneous Tracking (BEST) and involve simultaneous ranging from several 
sites. Although absolute station positions may be contaminated by unaveraged data 
and force model errors, these errors cancel in baseline calculations and therefore 
do not bias the rate-of-change estimates (8). 

The use of two different spacecraft results in an inhomogeous set of baseline 
estimates, which complicates the retrieval of an average rate-of-change for the 
entire time period covered by the data. Another source of inhomogeneity stems 
from changes in the locations of ground stations during the experiment. The 
southern site was changed from Otay Mountain near San Diego, California, to 
Monumeat Peak, ~50 km away, in 1981, and the northern site, located near Quincy, 
California, was moved by about 500 m in 1980. An adequate ground survey tying 
the two locations with an accuracy comparable to the SLR measurements is avail- 
able in the latter case, but not in the former. We are therefore dealing with a 
data set properly described as the juxtaposition of four subsets, each internally 
homogeneous but having unknown relative offsets in baseline length. These subsets 
are listed in Table 1, and the individual data points are plotted in Pig. 2. For 
clarity of graphical presentation, the baseline for each subset is arbitrarily offset 
by 0.4 m. 


We have analyzed the observations presented in Fig. 2 under the assumption 
that the baseline rate-of-change remained constant, for the entire 11-year period 
covered by the data. The time-dependent baseline length for the ith data subset 
(i-1 4) is therefore given by a modeling equation of the form 

B-(t) - B - (t-t Q ) + Bf, (1) 


Table 1: SLR data sets used to derive the rate-of-change of the Quincy-San 
Diego baseline. 


lj Set 

Iino. 

Satellite 

Baseline 

Solution 

Type 

No. 

Pts. 

r.m.s. 

(mm)+ 

Cumulative 

Importance 

SI 

BE-C 

Quincy-Otay Mt. 

Special 

4 

89 

0.26 

wm 

LAGEOS 

Quincy-Otay Mt. 

Annual 

3 

37 

0.58 

KB 

LAGEOS 

Quincy-Mon. Pk. 

Annual 

2 

<1 

0.06 

m 

LAGEOS 

Quincy-Mon. Pk. 

BEST 

14 

43 

0.10 


t Misfit with respect to constant-rate solution of Fig. 2. 
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Fig. 2: Variations in the Quincy-San Diego baseline length from SLR meas- 
urements between 1972 and 1982. Constant-rate solution was fitted simul- 
taneously to the four data subsets listed in Table 1 and shown here with 
arbitrary 0.4 m offsets. Error bars are standard deviations obtained from 
individual baseline estimates. 


where the intercept B? is the absolute baseline length at the fiducial time t , and 
B is a constant rate-of -change. 

Because of the offsets in the observed baseline lengths, the B^'s as well as B 
must be considered a priori unknown; hence, the data generate a system of 23 
equations in 5 unknowns. The four intercepts are only of ancillary interest, so it 
is not necessary to solve for their values directly. We therefore applied an 
appropriate projection operator to the linear system to annihilate the dependence 
of the data on these parameters and then solved for 6 by least-squares. The 
weights used in the least-squares criterion were'-jchosen to be proportional to the 
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standard errors depicted on Fig. 2. These formal errors of estimation are 
appropriate as relative weights within each homogeneous data subset, but there 
remains a subjective choice regarding the weighting of different subsets with 
respect to each other. In view of the earlier discussion, we down-weighted the 
BG-C data relative to the LAGEOS data by a factor of two in order to reflect the 
greater reliability of the latter. This has the practical effect of nearly equalizing 
the weights of the most reliable BE-C points and the least reliable LAGEOS points 
in the inversion. 

The least-squares formulation permits calculation of "data importances", addi- 
tive quantities which represent the information contributed by each datum or data 
subset to the solution (9). In the present case, since we are estimating the single 
scalar quantity 6, the data importances of the projected system sum to unity. The 
cumulative data importances for the four data subsets are listed in Table 1. 

The rate-of-change estimate obtained by the least-squares procedure is 

B “ -61.7 + 8.7 mm/yr, (2) 

where the uncertainty is the standard error of estimation. The r.m.s. residual for 
this solution is 52 mm, and the weighted r.m.s. residual is 29 mm. As seen on 
Fig. 2. this solution satisfies all but five observations at the one-sigma level, and 
all but one of the BEST estimates at the two-sigma level. 

The misfit to the 1972 BE-C measurement is no cause for concern, since it 
was the first collected for the SAFE baseline and is the most uncertain. Its indi- 
vidual importance is only 0.04; that is, it contributes but a small fraction of the 
information used to estimate B. In fact, completely removing this first point from 
the data set merely changes the answer to -59.9 + 8.9 mm/yr, which is not 
resolvable from the value quoted above. 

The absolute baseline lengths retrieved from the BE-C data and the LAGEOS 
data for Quincy-Otay (subsets 1 and 2 of Table 1) differ by 0.76 m. While such a 
discrepancy was not unexpected in view of the differences in data reduction tech- 
niques. its size indicates BE-C force-model errors that were somewhat larger than 
anticipated. However, our claim that the rate is not contaminated by the bias in 
absolute baseline length is supported by the consistency of the baseline changes 
between 1976 and 1979 for both spacecraft. 

As for the poorly fitted BEST points, it must be noted that this data reduc- 
tion technique, while minimizing the effects of dynamic errors on the baseline esti- 
mate, is quite susceptible to intermittent anomalies which occasionally occur in the 
laser data. Its reliance on a limited set of strictly simultaneous observations (in 
reality, manufactured laser "normal points" at two minute intervals) allows little 
averaging of the errore arising from anomalous data, causing an occasional spurious 
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baseline estimate. Fortunately, including the three outliers evident in Fig. 2 does 
not affect the slope estimate derived from the BEST subset. This subset 
comprises the majority of the data points but does not dominate the estimation of 
8; in fact, owing to its short time span, its cumulative data importance is only 
0 . 10 . 

By far the most important group of observations is the set of three annual 
LAGEOS solutions for the Quincy-Otay baseline between 1976 and 1982 (subset 2); 
its cumulative importance is nearly 0.6. These data cover approximately 5 years of 
observations, so they provide a good constraint on the slope. Subset 2 alone 
yields a rate of -53.6 + 11.4 mm/yr, again not significantly different from Eq. 2. 

Tapley et ah (7) produced an independent set of eleven estimates for the 
Quincy-Monument Pk. baseline from bimonthly arcs between October. 1981 and 
November, 1983. For this two year period, they deduce a rate-of-change of -64 + 

9 mm/yr. Their data reduction technique is different from ours, and they extend 
the coverage by a full year beyond the BEST solutions shown on Fig. 2. It is 
encouraging, therefore, that their value is in good agreement with Eq. 2. 

Assuming the Quincy site lies on a rigid North America plate (NOAM) and the 
Otay Mt./Monument Pk. sites on a rigid Pacific plate (PCFC), we can compute the 
baseline rate-of-change using the PCFC-NOAM instantaneous angular velocity vector 
derived from plate-tectonic data. The eleven-plate model RM2, which satisfies a 
globally distributed set of marine magnetic rates, transform-fault directions and 
slip-vector azimuths (10), yields B ” -52.9 + 2.7 mm/yr, consistent at the one- 
sigma level with the SLR estimate of Eq. 2. This removes the discrepancy noted 
by Smith et al. (2), who derived a preliminary SLR rate of -90 + 30 mm/yr from 
the 1972-1977 BE-C observations. As shown in Fig. 2, the BE-C and LAGEOS 
points prior to 1979, taken alone, do indeed suggest a higher value, which could 
indicate some temporal variation in the apparent rate, but the hypothesis that the 
deviations from a constant-rate model are random statistical fluctuations cannot be 
rejected at even a low confidence level (11). 
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Appendix 2 


CONSTRAINTS ON WESTERN I.J - S - DEFORMATION 
FROM SATELLITE LASER RANGING 


Abstract; The rate of shortening between San Diego and Quincy, California, 
measured by satellite laser ranging (SLR) is substantially greater than that implied 
by the rate of slip measured on the San Andreas fault. When combined with other 
kinematical constraints on western U.S. deformation, these comparisons limi t the 
extension rate across the Great Basin to be less than 9 mm/yr and imply signifi- 
cant deformation of the Pacific plate west of the San Andreas. Along the central 
California margin, the integrated rate of deformation is estimated to be 18 + 5 
mm/yr and entails compression perpendicular to the San Andreas, as well as right- 
lateral strike-slip motion parallel to it. Within the framework of our tectonic 
model, the SLR data limi t the integrated deformation rate across the Southern 
California Borderlands to be less than 7 mm/yr. 


The rate-of-change of the Qnincy-San Diego baseline (Fig. 1) measured in 
NASA’s San Andreas Fault Experiment (SAFE) since 1972. and more recently by - 
NASA's Crustal Dynamics Project is (1) 


6 - -61.7 + 8.7 mm/yr, (1) 

The value calculated from the plate-tectonic model RM2 (2) is -52.9 ± 2.7 mm/yr, 
which assumes Quincy lies on a rigid North America plate (NOAM) and San Diego 
on a rigid Pacific plate (PCFC). The geodetic estimate of the SAFE baseline 
shortening rate observed over an eleven-year period is thus compatible with plate 
motions averaged over a much longer interval (1-3 million years in the case of 
RM2). However, neither is easily reconciled with the hypothesis that PCFC-NOAM 
relative motion is localized on the San Andreas fault. . We consider two fiducial 
points juxtaposed on either side of the fault zone, r£ just to the west of the San 
Andreas and just to the east, and let Vg^ be the (steady-state) velocity vector 
of the former with respect to the latter. In central California, where the fault 
best approximates a simple transform, a range of geological and geodetic observa- 
tions local to the fault zone provides the estimate (3, 4) 


f SA 


f 34 + 3 mm/yr 

< 

l N41*W ± 2* 


( 2 ). 


The RM2 PCFC-NOAM velocity vector.-, calculated at the fiducial location 
r c - 36*N, 12Q.6*W (Fig. 1) is ' & - - ^ 
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original pmE>m 

OF POOR QUALITY 



Pis. 1: Oblique Mercator projection of the western United States about the 
RM2 PC PC- MO AM pole of rotation, showins major Quaternary faults and 
reference points used in oar vector calculations. C is the fidncial point on 
the San Andreas fault where the vectors in Bqs. 2-6 are specified. AB is 
the line alone which the expansion rate a£ the Great Basin is defined. 
Triangles are SLR stations at Quincy (Q) and San Diego (S), California. 
Pauits are from (5) and (6); lettered features are the Gariock fault (G), 
Hosgri fault (H),' San Andreas fault (SA), San Gregorio fault (SG), Sierra 
Nevada front (8N), Wasatch front (WF), and Walker Lane (WL). 
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r 55.9 ± 2.7 mm/yr 

Vp»i ” < (3) 

l N35.5*W t 1.9* 


The discrepancy between Eqs. 2 and 3 is significant at the two-sigma level in both 
rate and azimuth; it is thought to be a measure of the present-day tectonic defor- 
mation integrated across a broad zone of PCFC-NOAM interaction that extends from 
the continental margin of California to the Rocky Mountains (2. 5). ' 


How the deformation represented by this discrepancy is partitioned east and 
west of the San Andreas is a question of considerable geological interest and 
seismogenic implications. If Vq^ is the velocity of the fiducial point r£ in the 
reference frame rotating with NOAM and v^p is the velocity of rj~ in the PCFC 
frame, then the discrepancy vector Vpj| • v^, by definition, must equal 
Vqjj • v ( ..p: the deformation distributed over a plate boundary zone must sum up to 
the relative plate velocity, v^ thus summarizes the deformation occurring on the 
NOAM side of the San Andreas, whereas v^p summarizes the deformation on the 
PCFC side. We have formulated these vectors in terms of frame-independent 
integrals of the velocity gradient tensor along paths across the zone of deformation 
(3). 


The distribution of Quaternary faulting shown in Fig. 1, as well as the 
observed seismicity, defines an essentially rigid Sierra Nevada-Great Valley block 
(SNGV) extending from the San Andreas eastward to the Sierra Nevada front (5). 
The relative motion between this block and NOAM is taken up primarily by defor- 
mation in the Great Basin of Nevada and western Utah. Stress and strain orienta- 
tion data pertaining to Quaternary deformation of the Great Basin constrain the 
direction of this motion to be approximately N60*W (7). These constraints have 
been used to construct a one-parameter family of kinematical models parameterized 
by the total rate of expansion across the Great Basin from south-central Utah to 
northern California (along the line AB in Fig. 1){3>. We define a normalized velo- 
city vector v^jj by the equation 


v CN {«> - « v CN , (4) 

where oc is a dimensionless parameter such that ex X 10 mm/yr is the opening rate 
for this path. In other words, v^ is the unique velocity which would describe 
the local SNGV-NOAM motion at if the Great Basin were expanding at a total 
rate of exactly 10 mm/yr along AB. As documented in (3), the data compiled by 
Zoback and Zoback (7) yield the conditional estimate 
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a 

v 


CN 


f 10.1 + 0.7 mm/yr 
< 

[ N63*W + 5* 


( 5 ) 


The estimation procedure, which is based on a numerical search over possible 
SNGV-NOAM pole positions, explicitly accounts for sphericity of the earth's sur- 
face (on a flat earth, the modulus of Vqjj would be exactly 10 mm/yr). 

The locus of Vqjj parameterized by at is compared with the discrepancy vector 
Vppj - Vg^ in the upper vector diagram of Fig. 2. Clearly, the opening of the 
Basin and Range is not in the right direction to account entirely for the 
discrepancy, regardless of the value of ac. This implies that some deformation 
must be occurring west of the San Andreas fault. 

Given the definitions embodied in Eq. 4 and the constraints imposed by the 
plate-tectonic boundary conditions, the vector summarizing deformation west of the 
San Andreas can be recast in the form 


v cp (oc) “ « v CN - (v pN - v SA ), (6) 

and an estimate of v^p conditional an ac can be derived from Eqs. 2, 3 and 5. 
This expression and the corresponding equation for the variance matrix of Vgp(ac) 
yield the lower vector diagram in Fig. 2. The locus of v^pla) is a straight line in 
the direction of v^, and its two-sigma error ellipse sweeps a band approximately 
12 mm/yr wide, broadening slightly with increasing ac. For values of oc near zero, 
VQp(ac) is subparallel to the San Andreas and could be accommodated by predom- 
inantly right-lateral faulting localized along, say, the San Gregorio-Hosgri fault sys- 
tem, which parallels the California coastline (Fig. 2). On the other hand, if the 
Great Basin extension rate were as great as 20 mm/yr (ac - 2), as some students 
of Basin-and-Range deformation have advocated (9), then v^p would be essentially 
perpendicular to the San Andreas and would imply compression across the California 
margin with very little strike-slip motion. Interpretations of recently published 
seismic reflection profiles traversing the southern part of the Hosgri fault favor 
this model (10). 

The SLR observations discussed in (1) provide an upper bound on ac that 
excludes this latter possibility. Assuming the Quincy site lies on SNGV and the 
San Diego sites (Otay Mt. and Monument Pk.) on PCFC, we can derive an estimate 
of B conditional on a by the same procedure used to obtain Eq. 5. The result is 


6 - -(52.9 - 7.9a) + (7.29 + 0.49a 2 )* mm/yr. 


(7) 


- 16 - 



ORIGINAL PACE" IS 
OF POOR QUALITY 



PIC. 2z May Lambert conformal conic projection of central California show* 
ing major zones of Quaternary faulting (8). Lettered fault zones are the 
Big Pine (BP), Calaveras (C), Hosgri (H), Nachaiento (N), Pilarcitos (P). Rin- 
conada (R), San Andreas (SA), San Gregorio (SG). Santa Lucia Bank (SLB), 
and Santa Ynez (SY). C is the fiducial point at 36*N on the San Andreas 
used in vector calculations. Upp e r inset: Velocity diagram in the tangent 
plane at C showing discrepancy vector Vp^ ■ and vector Tqj^oc) 

describing Basin- and- Range defonsatioa; a - 1 corresponds to an opening 
rate of exactly 10 mm/yr across the Great Basin along line AB in Fig. 1. 
Error ellipses are two-sigma. Lower Inset: Velocity diagram at C showing 
?0p ( which describes defonsatioa west of the San Andreas. Dathed line is 
locus parameterized by a. Error ellipses are two-sigma. 


- 17 - 


Western U.S. Deformation from SLR 


As ac increases from zero, the total expansion of the Great Basin accelerates, and 
the agreement between Eqs. 1 and 7 degrades. If ac were as great as two, for 
example, then the value of the SAFE baseline rate-of -change would be -37.1 + 3.1 
nun/yr, which is not consistent with the SLR results. 

A formal upper bound on at can be obtained from Eqs. 1 and 7 by applying a 
one-sided chi-square test; at the 95% confidence level, we find at < 0.80. If 
instead of Eq. 2 we use the estimate of -64 + 9 mm/yr given by Tapley et al, 
(11), this bound is reduced to 0.56. The chi-square test assumes the standard 
deviations assigned a priori to the SLR data are exact; scaling these uncertainties 
to the root-mean-square misfit provided by the constant-slope solution reduces the 
standard error of 6 in Eq. 1 to 8.3 mm/yr, and applying the appropriate one-sided 
F test gives ot < 0.89. 

The SLR measurements thus limit the extension rate of the Great Basin to be 
less than 9 mm/yr. This is consistent with, but more stringent than, the upper 
bounds obtained from direct geological and geophysical observations of Basin-and- 
Range deformation (12). Moreover, through Eq. >7, the SLR data constrain the 
vector V£p to have a significant projection parallel to the San Andreas; the lower 
bound on this component of motion is 5 mm/yr, and it is likely to be considerably 
larger (Fig. 2). In fact, our results are consistent with strike-slip rates of 6-17 
mm/yr inferred from post-Pliocene offsets along the San Gregorio-Hosgri fault sys- 
tem (13). Hence, a model with purely compressive deformation along the California 
margin is excluded. We note, however, that ac is constrained by observations in 
the Great Basin to be greater than 0.1 (12), requiring the component of v^p per- 
pendicular to the San Andreas to be greater than zero, as implied by geological 
observations of regional compression west of the fault (10). 

We infer that the rate of deformation west of the San Andreas is large. For 
ex - 0.9, our best estimate of v^-p has a magnitude of 16 mm/yr, and it is 
increases to 21 mm/yr for oc - 0.1. These estimates are one measure of the 
3eismogenic potential of the central California margin. 

We have assumed throughout this analysis that the San Diego sites are rigidly 
attached to PCFC. Any geologically plausible deformation west of these sites -- 
accommodated, say, on offshore faults within the seismically active California 
Borderlands- would imply a decrease in the magnitude of B relative to Eq. 7 and. 
hence, would lower our upper bound on ac. For example, the block-tectonic model 
of southern California developed by Bird and Rosenstock (14) predicts that the San 
Diego sites move with respect to PCFC at ~3 mm/yr in the direction ~N40'W; 
such a motion would lower our bound on ac from 0.9 to about 0.6. Conversely, a 
lower bound on at can be used to place an upper bound on the rate at which the 
San Diego sites translate with respect to PCFC. Assuming this motion is N40*W, 
and adopting at > 0.1, we find that this rate must be less than 7 mm/yr to be 
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consistent with the SLR data. The integrated rate of deformation across the 
Southern California Borderlands is therefore inferred to b<* less than that across 
the central California margin. Presumably, the difference between these deforma- 
tion vectors is taken up by the northward compression of the western Tranverse 
Ranges (14). 

The preliminary results in this report point the way to future applications of 
SLR and other space-geodetic techniques in the study of regional deformation. In 
particular, they illustrate the advantages of combining SLR measurements with 
other geological, geophysical and geodetic data sets. When considered in the con- 
text of regional kinematical models constructed from such data sets, the SLR 
observations from a single California baseline yield quantitative constraints on the 
distribution of deformation in the western United States. Much more stringent 
constraints can be expected from the large set of North America baselines 
currently being measured under the auspices of NASA's Crustal Dynamics Project 
(15). Space-geodetic observations should be especially useful in fixing the rates of 
regional deformation that are difficult to measure with conventional techniques but 
are critical to the assessment of seismogenic potential. . 
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Abstract 


Satellite laser ranging observations to Lageos during 1978-1983 are used to 

analyze the vector rate-of-change of annual station positions for comparison with 

global plate tectonic models. We develop an algorithm suitable for the reduction of 

large, inhomogeneous space-geodetic data sets, which accounts for the full position 
covariances, and coordinate transformations between epochs. Tectonic motions are 
isolated by reducing systematic errors in radial position determination and 
contaminations due to errors in earth-orientation parameters. This is accomplished 

by forming a linear system of equations which describe differenced station position 
vectors as a combination of tectonic motion model parameters, and nuisance 
parameters which comprise non-tectonic motions such as rigid-body rotations. The 
system is normalized with the square-root of the covariance matrix, and then 

rendered insensitive to the nuisance parameters by the application of projection 

operators. Model parameters are estimated from the resulting reduced system by 

generalized-inverse techniques. This methodology is applied to the problem of rate 
scaling of global plate tectonic models. The model is parameterized by ym 0 , where m Q 
is an a priori model of tangential relative velocity vectors, such as RM2; thus, y = 1 
implies that the mean rate of motion of the observing stations is identically equal 

to the RM2-predicted rates. We find y = 0.95 ± 0.16, indicating that the mean rate of 

SLR-derived motions for a five-year interval are consistent with RM2, which is 

averaged over a 2-m.y. interval. We find that the data do not place significant 

constraints on possible errors in marine magnetic anomaly timescales, but are not 

consistent with an interval of quiescent or accelerated plate motions. 
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Introduction 


Satellite laser ranging (SLR) and very-long baseline interferometry (VLBI) are 

space-geodetic techniques currently capable of determining relative locations of 

widely separated positions on the earth's surface to within sub-decimeter accuracy 
(Shapiro, 1983). Results from interferometry with the Global Positioning System 
(GPS) are achieving similar accuracies (Bock et al., 1985). This capability raises the 
possibility of detecting and monitoring relative plate motions and deformations. 
Analysis of the rate of change of baseline length between many station positions are 

in good agreement with predictions based on global plate tectonic models (e.g. Smith 

et al., 1979; Christodoulidis et al., 1985; Tapley et al., 1985; Clark et al., 1985), although 


some significant discrepancies 

do exist (e.g. 

Minster 

and 

Jordan, 

1984). 

These 

treatments of 

space-geodetic data concentrate 

on the 

scalar 

baseline 

lengths, 

which 

are relatively 

insensitive to 

systematic errors due 

to coordinate 

transformations 


between epochs. However, the data are intrinsically three-dimensional. Thus, if the 
errors can be removed or reduced in a systematic manner, an analysis of the rate of 
change of the geocentric position vector could provide additional insight, especially 

with regard to directions of motion. As these data sets grow, the need for efficient 

algorithms capable of handling large, inhomogeneous, multiple-epoch data sets 
becomes increasingly critical. In this paper, we present the nucleus of such an 

algorithm, and then, using data derived from laser ranging to the Lageos satellite, we 
apply this algorithm to the problem of rate scaling of global plate tectonic models. 

We first examine the nature of the satellite laser ranging data and some of its 

inherent difficulties. In May 1976, Lageos (Laser Geodynamics Satellite) was 
launched into a nearly circular orbit at an elevation of 5900 km and an inclination of 

109.8° to the equatorial plane. The satellite is a sphere 60 cm in diameter, weighing 
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411 kg, and is covered with 426 retroreflectors. The satellite’s design and orbit 
minimize the fore'** due to solar radiation pressure, earth albedo, and atmospheric 
drag, and thus, Lageos very nearly defines an artifical reference frame in earth 

orbit. Ground-based laser systems measure the round-trip pulse propagation time 

and, using an adopted value for the speed of light, calculate the range. In theory, if 

the satellite's position is known perfectly, the location of the laser site may be 
determined after a number of ranges have been obtained. In practice, of course, the 
problem is considerably more difficult. Since ranges to Lageos, in general, are not 
made simultaneously, the orientation of the earth with respect to the satellite 

reference frame must be determined for the time of the observation. This requires 
knowledge of a complex set of motions which include earth precession, nutation, 
polar motion, and rotation, as measured by UT1. The satellite's position is controlled 
by the external forces which act upon it; these include the earth's gravity field, 

luni-solar and planetary gravity fields, solid-earth and oceanic tides, solar radiation, 
and earth albedo. Range measurements are affected by refraction of the laser pulse 
in the atmosphere, tidal displacements, pulse detection accuracy, and are scaled by 
the adopted value for the speed of light. Thus, the accuracy of the station positions 

depends critically upon the accuracy of the earth-orientation, force, and 

measurement models. One such set of models and the procedure for reducing raw 

range values to station position coordinates used by Goddard Space Flight Center is 

called SL5.1AP, and is described in Christodoulidis et al. (1985), and Smith et al. ( 1985). 
The data that we utilize is derived from an improved set of models, denoted by SL6, and 
are solutions for annual geocentric vector positions and covariances derived from 
observations of Lageos between 1978 and 1983. 
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Methodolog y 


We now consider the problem of constraining plate tectonic motions and 

deformations from space-geodetic networks. Ideally, the analysis procedures should 
be capable of handling data derived from disparate techniques such as SLR, VLBI, and 
GPS simultaneously, taking into account the coordinate transformations and scaling 
parameters particular to each. For example, the natural coordinate system for SLR 
has an origin at the earth center of mass, whereas VLBI uses an origin at the solar 
system barycenter; also, SLR, unlike VLBI, is sensitive to the adopted value of the 
earth gravitational mass. Any estimation of parameters by inversion techniques 

should account for possible errors in the data, therefore the analysis must 

incorporate the full covariance matrix corresponding to the data. Another potential 
problem arises from non-uniform sampling of network stations from epoch to epoch. 
Usage of homogeneous data sets, i.e. including only stations that are sampled during 
every epoch, is often too restrictive and ignores important sources of information. 
However, inhomogeneous sampling may introduce biases in the estimation of station 

parameters due to changing network geometry, which result in unequal weighting 
between successive epochs. For example, one epoch may preferentially sample "fast" 

plates. The approach we have taken is similar to that used by Jordan and Sverdrup 

(1981) for the estimation of earthquake location, whereby we form a linear system of 


equations 

which 

describe 

differenced 

station 

position 

vectors as 

a 

combination 

of 

tectonic 

motion 


model 

parameters 

and 

nuisance 

parameters, 


which comprise 

non-tectonic motions 

such 

as rigid-body 

rotations. 

The system 

is 

normalized 

with 

the square-root 

of 

the 

covariance matrix. 

and then 

rendered 

insensitive to 

the 

nuisance 

parameters 

via 

application of 

projection operators, a process 

which we 

term 


"denuisancing” . 
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Let us consider a simplified algorithm tailored to the specific problems presented 
by SLR observations. First, if all the data is processed in a consistent manner, such as 
with the SL6 model, scaling of the network due to the adopted value for the speed of 
light will be uniform and hence negligible. Similarily, the earth center of mass, 
which is an implicit function of the assumed gravity model, is constant for all epochs. 
Thus, representation of scaling or transformation of coordinate systems may be 

safely ignored. Errors in earth-orientation parameters, such as polar motion and 
UT1, will result in rigid-body rotation of the entire network about the earth center of 
mass. However, absolute plate motions, which are relative plate motions defined in a 
particular angular velocity reference frame, will tradeoff exactly with rigid-body 
rotations. Therefore, we may only resolve relative plate motions from the data and 

will treat rigid-body rotations as nuisance parameters. Some provision should be 

made for differences between errors in site motion tangential and radial to the 

surface of the earth, which arise from inhomogeneous sampling of ranges. For 

example, lasers may range at all horizontal azimuths, but may not range at vertical 
angles below the horizon; thus, systematic errors in satellite location and in the 

atmospheric refraction model are expected to map primarily into the radial 
components of position. Therefore, radial components may also need to be treated as 
nuisance parameters. Finally, tectonic motions may bias the denuisancing of 

rotations, since incomplete sampling of the plate system might introduce rotations of 

tectonic origin. This subtle problem may be minimized if prior information about 

the motions is known; the a priori tectonic model is subtracted from the data before 
denuisancing, and then is added to the resulting model after the inversion has been 
completed. This process, which we term "bias stripping", causes the component of the 
solution which is not constrained by the data to coincide with the a priori model, thus 
care must be exercised in the choice of the a priori model. 

With these considerations in mind, we may formulate a specific approach. We 
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define x to be a vector composed of the three-component geocentric vector station 

positions for all stations sampled during all epochs under study. These positions will 

A 

have associated covariances defined by V x = <(x-<x>)'(x-<x>)^>, where the brackets 

A 

denote expected value, and the T transpose. x and V x are the quantities supplied by 
the SL6 solution. This covariance matrix is not an accurate reflection of the true 
errors associated with the station coordinates due to systematic errors in the normal 

point ranges used to derive the SL6 solution. Normal points are formed by grouping 
raw ranges to Lageos from a laser site into 2 minute bins, and then finding average 
ranges from those samples; this data compression, which can typically be two orders 
of magnitude, is useful computationally and reduces biases due to uneven sampling of 
the orbit of Lageos. The precision of ranges is on the order of 2-9 cm, and thus the 

formal precision of the normal points is generally on the order of 1 cm or less. 
However, due to systematic orbit model errors, which propagate as the square of time, 
the actual accuracy of the normal points is more on the order of 10-15 cm over a 30 
day orbital arc. To reconcile this problem, SL6 uniformly assigns each normal point a 
formal uncertainty of 100 cm, which means that the stated covariances are probably 

an overestimation of the errors. A more realistic estimate of covariance might be 

A 

V x = v 2 V x , where v 2 is a scalar to be determined a posteriori by the 
sampling variance, which is a function of the weighted error vector and the number 

of degress of freedom. 

We define 

d = E - x (la) 



where E is the operator composed of l’s, 0's, and -l's which difference vector 
positions between successive epochs for all sampled stations. d therefore represents 
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Based on the 


the apparent vector motions of the stations from epoch to epoch, 
previous discussion, we consider a decomposition of these motions into rigid-body 
rotations, tectonic motions, such as rigid plate motions and internal plate 

deformations, radial motions and errors. We construct a linearized set of equations of 
the form: 


A ' m + B • n = d (2) 

where m is a vector composed of model parameters corresponding to tectonic 
motions, and n is a vector composed of nuisance parameters. The particular 

representation of these vectors will depend on the desired parameterization. For 

example, if rigid-body rotations are to be denuisanced from the system, n will be 
composed of the three-component angular velocity vectors corresponding to the 
rotations between successive epochs, and B will be composed of the 3x3 
anti-symmetric infinitesimal rotation matrices times the epoch-time differences for 

each station sampled. m may be the denuisanced vectors of motion, in which case A 
will be the identity matrix; or, as discussed in the following section, m may be a 

rate-scaling parameter and A the vectors corresponding to an a priori global plate 
tectonic model. Explicit representations of these matrices will be given for the 

specific application of rate scaling of global tectonic models; however, we wish to 
emphasize that this is a completely general formulism capable of handling many 

possible parameterizations, provided the system is linear or can be linearized. 

We normalize the system by weighting with the square-root of the covariance 

matrix 

A “ A A 

A • m + B • n = d (3) 

A A 

where A= V d ' 1/2, A, etc.; note that the carat now indicates normalization by 
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1/0 

V d . The system is denuisanced by applying 

A A A 

Q B = I - B* B', (4) 

where the dagger, denotes the Moore-Penrose generalized inverse (Penrose, 1955). Q B 

A 

is the projection operator which projects (3) into the null space of B T , thus 

rendering the system insensitive to the nuisance parameters; it is both symmetric 
and idempotent (Q B = Q B T = Q B ‘Q B )• Multiplying both sides of (3) by (4), and making 

A 

use of the fact that Q B • B = 0 , we find 


A A 



Q B * A- m = Q b - d. 

(5) 

The 

model parameters are estimated by taking the generalized 

inverse of the 

hand side 

of the reduced system (5) 

* A x A 



m = (Qb" a) Q b * ^ > 

(6) 

with the 

associated unsealed covariance matrix 

A A x A ^ 



V m = a ) ,, Qb , ^B* A ^' 

A A A 

(7a) 


= ( A • Q g • A ) . 

(7b) 


7 

We now estimate the covariance scaling parameter v ~ by comparing the misfit 
between the model and the data; this is defined by the normalized error vector, which 
is a random variable found by substituting (6) into (5): 


A - A A 


= Q B - d - Q b " A m 

(8a) 

A A x A 

= (I ■ (Q g ' A )’(Q b ' A) )• Qg- d. 

(8b) 
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The squared length of the normalized error vector divided by the covariance scaling 


parameter is itself a random variable and is chi-squared distributed with N-M degrees 


of freedom, where N is the rank of the reduced data space Q Q * d , and M is the rank of 


A A , 

the projection operator (Q B ‘A) , (Qg' A)' (e.g., Jordan and Sverdrup, 1981). Thus, it 
has the expected value 


< e T * e> / v 2 = N-M. (9) 


The observed misfit is one realization of this error process, and can thus be used to 
estimate v 2 . This estimate is called the normalized sampling variance: 




/(N-M). 


( 10 ) 


Application of bias stripping is a simple extension of the method outlined above. 
Let m 0 be the a priori model which is assumed to be known perfectly; that is, let m 0 
have zero variance (m Q = < m Q >). The a priori model is removed from the linear 

system (2): 


A'(m-m 0 ) + B*n = d - A • m Q . (11) 

The equations are normalized with the square-root of the covariance matrix in (1), 
since m 0 has zero variance, and applying the appropriate denuisancing and 
generalized-inverse techniques, we find the estimate of the model parameters 

^ A 4. A A 

m = m Q + (Q B ’ A) ' • Q B ' ( d - A - m 0 ) (12a) 

= (Qg- A)' - Qg- d + (I - (Qg-A)'-(Qg - A))-m 0 . (12b) 
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It is straightforward to show that the unsealed covariance matrix, the normalized 
error vector, and the sampling variance are given by (7), (8), and (10) respectively; 
this follows from the fact that m 0 has zero variance. The first term in (12b) is 
equivalent to (6), the result obtained with bias stripping. The second term projects 
m 0 into the model space which is unconstrained by the data. Thus, the component of 

the model which is not constrained by the data is made to coincide with the a priori 

model. 

We have described a very general approach for the estimation of model 
parameters and variance scaling in the presence of nuisance parameters. We now 

apply the methodogy outlined above to the specific problem of the rate scaling of 

global plate tectonic models. 

Rate Scaling of Global Plate Tectonic Models 

As an introduction to this problem, we briefly describe one successful model of 

plate tectonic motions obtained by Minster and Jordan (1978), hereafter denoted by 
RM2. RM2 is a relative motion model based on the assumption that the earth's surface 
is composed of eleven plates which are rigid spherical caps moving tangentially to 

the surface. Directions of motion are estimated from azimuths of transform faults and 

earthquake slip vectors, and rates of motion are estimated from marine magnetic 
anomalies. No data are used from locations where the assumptions do not appear to be 

valid, such as regions of deformation at continent-continent boundaries. Magnetic 

anomalies 2, 2', and, occasionally, 3 are employed, thus the obtained rates are 

averaged over a 2-m.y. interval. One-sigma errors for the rates are on the order of 5 
to 10%. 

As an application of the algorithm decribed in the previous section, we estimate 
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from the SL6 solution the RM2 rate-scaling parameter y. That is, we assume that the 

orientation of the relative angular velocity vector for each pair of plates is known 

perfectly and that the rate of rotation is of the form yco, where to is the RM2 angular 
rate of rotation. Thus, the relative magnitude of the RM2 rates are assumed to be 

known perfectly, but they may be scaled according to y. For example, y = 1 implies that 
the relative velocities of the observing stations are identically equal to the 

RM2-predicted values, whereas y < 1 implies that the observed mean rates are slower 

than RM2. Thus, we seek to characterize the data by a single model parameter. Such 
a characterization maximizes the resolution of the model parameter at the expense of 
greater variance. As such, it represents a logical first assessment of the data, and will 
help to ascertain the appropriateness of RM2 as an a priori bias-stripping model. The 

rate-scaling parameter y is a useful description of mean global plate motion, and 
facilitates comparison of mean rate estimates derived from geophysical and geodetic 
techniques representing timescales over many orders of magnitude. The utility of 

such comparisons will depend on the assumption that the relative directions of plate 
motion do not change substantially over the averaging interval. For example, 

estimations based on plate reconstructions over the last 80 m.y. (e.g., Davis and 
Solomon, 1981) which include major reorganization of plate geometry, are possible, 
but will be degraded due to poorly-matched motions. RM2 provides a reference 

estimate of y based on marine magnetic anomalies averaged over a 2-m.y. interval, 
whereas space-geodetic data offer independent estimates which sample the last 10 

years. 

We now adapt the methodology of the previous section to the problem of 

estimating y from SLR data. Let d be the vector composed of the SL6-derived 

• A 

geocentric vector station motions, with associated covariance V d , as given by (1). To 
compare d to RM2-predicted tangential motions, we must reduce biases due to 

rigid-body rotations and radial motions from the data. An infinitesimal geocentric 
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rigid-body rotation may be represented by p*£5, where £1 is a three-component 
rotation vector, and p is the 3x3 anti-symmetric matrix 


P = 


0 z -y 

-z 0 x 

y -x 0 


( 13 ) 


where x,y,z are the nominal cartesian station coordinates (e.g., Goldstein, 1980). We 

define B to be composed of the p ; ’s corresponding to motion of the ith station, and n 
be composed of the flj k 's corresponding to motion between the jth and kth epochs; in 
this way, the rigid-body rotations described by n may be treated as nuisance 
parameters. The assumption of infinitesimal rotation is valid for the observed 

centimeter motions; such rotations are additive, allowing data from inhomogeneously 
sampled networks to be fully utilized in a simultaneous inversion over all sampled 
epochs. For example, for a two station network where station 1 is observed during 
three epochs, and station 2 is observed in the first and last epochs, we would define 


B • n = 


Pi 0 

0 Pi 

P2 P2 



In 


23 J 


(14) 


To remove radial components of motion and all covariances between tangential 

and radial motions, we convert from geocentric to local geodetic coordinates via the 

transformation matrix J, which has jacobian block diagonal elements corresponding 
to the nominal station positions, and isolate the local north and east components of 

motion by removing the rows and columns corresponding to vertical components via 

the tangential selection matrix T, which is simply composed of l's and 0's. 

Thus, we define 


d' = F- d 


(15a) 
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(15b) 


V d , = F* V d -F T 
B' = F- B (15c) 


where F = T* J. Because of the linear relationship between local geodetic velocities 
and angular rate of rotation (Minster et al., 1974), we may parameterize the model by 
ym 0 , where m Q is the vector composed of predicted RM2-derived velocities in local 
geodetic coordinates multiplied by the time difference between the epochs. The RM2 
velocities are calculated with respect to a fixed plate; the choice of the fixed plate is 
arbitrary since it results in a rigid-body rotation the network which is removed by 
application of Q B , as given in (5). Thus, we may define the system: 

Y m Q + B'- n = d'. (16) 

Noting the equivalence of (16) to (2), where y = m , and m 0 = A, we find from (6): 

A A 

y = (Q B '- ,n o) t ’ Qb' - d ' U?) 

which, being scalar, has unsealed variance, given by (7) 

CT y 2 = ( m o T 'QB' m o )t < 18 > 

and v 2 is estimated from (8) and (10). 

This method does not account for possible biases in model estimation due to 
tectonic motion; we previously suggested that such contaminations might be reduced 
by bias stripping with an a priori model of tectonic motion. It is straightforward to 
show, however, that this procedure is applied implicitly when estimating the 
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rate-scaling parameter. Define y(a) = 8y + a, where a is a constant. We bias strip using 
am 0 as the a priori model; thus, a = 0 corresponds to no bias stripping, and a = 1 is 
equivalent to* bias stripping with m Q . Then (16) becomes 


It follows that 
and 


8y m Q + B'- n = d' - a m Q (19) 

8y=(Q B ,-m 0 ) r - Q B ,- d' - a (Qg,* m 0 ) t -Qg,*m 0 (20) 

y(a) = (Qg,* m 0 ) f - Qg,* d' + a [1 - (Qg,- m Q ) t - Qg,- m Q ] (21) 


Note that (Q 0 ,* m 0 )^* Qg,* m Q * s an idempotent scalar projection operator, and must 

A _ 

equal either 0 or 1. Qg,* m Q is the relative motion model in a no net-rotation frame; 
for non-trivial cases this will be non-zero and therefore the projection operator 

equals 1. The quantity in the square brackets of (21) vanishes and we find the 

estimate of the rate-scaling parameter is not a function of a; hence, bias stripping is 

implicitly contained in this particular formulation. 

The nominal station positions (shown in Figure 1) and distribution of normal 

points for the stations used in the estimation of the rate-scaling parameter are listed 
in Table 1; the inhomogeneous sampling of network geometry which is evident from 

this distribution provides some justification for the multiple-epoch denuisancing 

procedures described above. The data, which consist of 52 local geodetic motion 
vectors corresponding to 23 stations sampled over 6 epochs (1978-1983), are listed in 
Table 2; coverage includes 5 of the 11 RN12-defined plates. These vectors with their 

associated tangential error ellipses (1 sigma, unsealed) are shown in Figure 2; this 

figure provides some sense of the distribution and quality of the formal errors of the 
data--it does not portray the often significant covariance between station motions. 

Applying the above methodology to the complete five-year interval data set, we 
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obtain: 


y = 0.95 ± 0.16 (one sigma). 


The one-sigma level of confidence is obtained from s 2 = 0.355, and cr^ 2 = 0.072 . 
Additional estimations based on epoch-subsets are listed in Table 3; for example, using 
all the stations with data in both 1982 and 1983, which we denote {5,6}, the one-year 
interval estimate of y is 0.51 ± 0.50. These results are plotted as a function of the 
number of data vectors utilized in each solution in Figure 3. In all but four cases, the 
SLR-derived rate-scaling parameter is consistent with the RM2-reference value (y = 
1) at the one-sigma level of confidence. One exception, {3,4,5}, is very nearly 
consistent (1.275 ± 0.270), and the remaining exceptions ({1,2}, {1,2,3}, {1,2,3, 4}) all 
include data from the first epoch, 1978. This epoch was sparsely sampled, using 
measurement techniques which have been subsequently refined; it is therefore 
possible that these data may be subject to systematic errors not accounted for by SL6. 
If the 1978 data is ignored, estimates of the rate-scaling parameter averaged over a 


one-year 

interval 

vary 

between 

0.44 

and 1.31, and 

show 

less 

variability 

for longer 

intervals. 

In summary, 

we find 

that 

the SLR-derived 

mean 

rate 

of motion 

among the 

observing 

stations 

for 

a five-year averaging interval. 

and 

even 

for intervals as short 


as one year, are consistent with RM2. 

Discussion 

This result supports previous scalar baseline rate-of-change analyses indicating 
that the direct measurement of plate tectonic motions is feasible using space-geodetic 
techniques (e.g., Christodoulidis et al„ 1985); however, we find that an analysis of the 
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vector motions associated with network station position provides additional 

information and constraints on "instantaneous" plate motions. For example, the use 
of RM2 as an a priori model of tectonic motion appears to be reasonably justifiable, as 

the improvement in results for longer averaging interval indicates. It should be 
noted that estimation of the rate-scaling parameter assumed that RM2 is known 
perfectly with the exception of the absolute magnitude of the rates. In reality, RM2 is 

a model with associated covariance due to errors in azimuths and magnetic anomalies; 
however, the errors associated with RM2 are small compared to those of SLR data, and 
therefore are not likely to make any substantial contribution to the estimation of the 
rate-scaling parameter. Further experimentation with alternative a priori models 
which may provide additional improvements, such as NUVEL-1 (Demets et at., 1985), 
should, of course, be pursued. As data coverage becomes more extensive, global plate 

tectonic models based solely on space-geodetic data may ultimately supplant those 

based on less instantaneous data; these improved models will help to reduce biases in 

the estimation of model parameters. 

Another notable feature of this vector analysis is the estimation of the sampling 

A 

variance s 2 = 0.355. Previous estimates of the sampling variance, which have been 

on the order of 0.10 to 0.15 (D. Smith, private communication), have relied on the 

analysis of single scalar baseline lengths, under the assumption that additional 

variances corresponding to other baselines, even if they are larger, will not increase 

the overall variance of the system. This assumption is valid provided the variances 

are uncorrelated. This is not the case for the SLR data; use of the full covariance 

matrix as well as a unified algorithm which considers all the data simultaneously 

indicate that the previous estimates of sampling variance are too small by a factor of 

2 or 3. 

We now consider some of the broader geophysical implications of these results. 

We have found that the SLR-derived mean global plate rates for a 5-year interval 
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agree with the 2-m.y. average RM2 rates; this suggests possible limitations for the 

temporal variability of mean global motion. Anderson (1975) proposed that observed 

tectonic motions in the vicinity of plate boundaries, such as subduction zone 
trenches, may be relatively stationary for long periods (30-100 yr), except for brief 
(5-10 yr) "breaking cycles" in which the accumulated stress along the boundary is 
released. During such a cycle, the plate motions near the boundary are accelerated. 
Episodic motion at plate boundaries is, of course, well documented; Anderson, from 
observations of seismic migration patterns such as the Aleutian sequence between 

1899 and 1905, suggested these episodic motions are correlated and that the 
compression of the oceanic lithospheric plate between breaking cycles may extend 
well into the interior of the plate. One possible implication of this model is that global 
plate motions may be "jerky”; that is, quiescent except for relatively short periods of 
accelerated motions. The five-year SLR-derived mean rate of motion is not consistent 
with either a quiescent, or accelerated interval of plate motions. One-year estimates 
display graeater variability but are, with the exception of the 1978 data, also 

consistent with RM2. Possible interpretations of systematic temporal variations must 
be regarded with some skeptism due to the sparcity of data, and to the 
previously-mentioned systematic improvements in data collection techniques which 
may not be accurately modeled by SL6. Clearly, many additional years of space- 
geodetic data collection are required before strong constraints may be placed on 
temporal variation in mean plate motion. 

The sampling frequency of the SLR data severely restricts efforts to make 

interpretations regarding longer averaging intervals; this point is illustrated in 

Figure 4 which shows estimates of global rate-scaling parameter as a function of the 
sampling interval. Space-geodetic data provide the most recent estimates averaged 
over the shortest intervals; it is easily possible to imagine temporal variations of 
motion over the entire spectrum of timescales out to 2 m.y. or more. For example, 
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Vogt (1985) suggests, from a detailed examination of the magnetic lineations at six 
accreting plate boundaries, that opening rates have been up to 7-35% faster over the 

last 1 m.y. than over the period 1-2 m.y. Other geological techniques may provide 

crude estimates _ for timescales intermediate between space-geodetic and marine 
magnetic anomaly intervals. Estimates of slip rate at plate boundaries based on 

paleoseismicity, such as those obtained by trenching on the San Andreas Fault by 
Sieh and Jahns (1984), sample averaging intervals of 1000-10000 years. A robust 
estimate of y would require many such estimates of slip rate over the same time 
period; the natural paucity of suitable sites will hinder such an effort. Furthermore, 
the rates obtained by Sieh and Jahns (1984) are consistent with more recent geodetic 

estimates but show a large discrepancy with RM2 rates; Minster and Jordan (1984) 

suggest that distributed deformation both east and west of the San Andreas may 
account for this discrepancy. Thus, care must be exercised when measurements 

obtained at plate boundaries are used to estimate global plate motion. Despite such 
obstacles, the variability of the rate-scaling parameter over different timescales 

warrants further investigation, because such variations might be a reflection of the 

dynamical processes which drive the observed plate motions. 

As has been noted, RM2 rates are derived from observations of marine magnetic 
anomalies; errors in the magnetic anomaly timescale will therefore result in scaling 

errors of the RM2 rates. Under the assumption that global plate motion has been 

nearly uniform over the last 3 m.y., space-geodetic data offer an independent 

constraint on possible errors in the magnetic anomaly timescale. RM2 uses the 

magnetic anomaly timescale of Talwani et al. (1971) for the estimation of the age of 

anomalies 2, 2’ and 3. In a more recent compilation of anomaly timescales, Ness et al. 

(1980) point out several disadvantages of the Talwani scale. First, it was obtained from 
profiles over the Reykjanes Ridge, which is a slow spreading ridge; such profiles tend 
to have poorer resolution due to contamination ot normal and reversed blocks from 
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subsequent flows. Second, it is a hybrid scale, constructed by assuming the age of the 
end of anomaly 5 (9.94 m.y.) of Heirtzler et al. (1968) to be correct, assigning residuals 
from the Heirtzler scale to be due to errors in that scale, and then compromising on 

the age of the 3.1' anomaly from considerations of magnetostratigraphy. Finally, it 
does not incorporate the change in K-Ar decay and abundance constants suggested by 

the International Union of Geological Sciences in 1977; this change alone would 
result in a 2.68% decrease in RM2 rates (Mankinen and Dalrymple, 1979). Ness et al. 
(1980) propose a new anomaly timescale which displays systematic differences with 
the Talwani scale. Table 4 compares the chronologies, and the corresponding 

rate-scaling parameter obtained by dividing the age from the Talwani scale by teh 

age from the Ness scale. These estimates of y are shown in Figure 4. The variability of 
these estimates provides some measure of the possible errors in the RM2 marine 
magnetic anomaly timescale; we find that the variations lie well within the formal 
error bars of the SLR-derived estimate. Although there are indications that the 
anomaly timescale used to determine the RM2 rates may be in error, the quality of the 
SL6 data set is not at present sufficient to place useful constraints on that error. 


We have developed a general methodology for the analysis of space-geodetic data 
with the goal of constraining plate tectonic motions. The generalized-inverse 
analysis of the vector rate-of-change of the geodetic network accounts for the full 
covariance of the station motions and reduces contaminations due to rigid-body 
rotations and poorly-constrained vertical motions. We have demonstrated the 

applicability of these algorithms to obtain the RM2 rate-scaling parameter, which we 
believe to be a useful description of mean global plate motion worthy of further 
investigation via other geophysical and geological techniques. An estimate of the 
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rate-scaling parameter from a five-year interval of laser ranging to the Lageos 

satellite is found to be near unity and consistent with RM2 rates. There is a 
suggestion that the magnetic anomaly timescale on which RM2 is based predicts 
faster rates than more recent timescales, but SLR data do not as yet place any 
constraints on this error. In addition, the mean rate of motion of the plates argues 

against the suggestion that the plates are in a period of quiescent, or accelerated 

motion. 

An obvious extension of this analysis would be the inclusion of VLBI data which 

presently offer useful network stations in Europe. More complicated 

parameterizations of the model are under consideration. Examination of vertical 

motions would provide a further test of consistency, and might indicate areas of 

active deformation. Grouping the stations by assumed plate might yield constraints 

on the directions of tectonic motions, as well as the stability of the plates themselves. 
Allowing individual stations to move with constant velocity might reveal anomalous 
behavior at plate boundaries. Space-geodetic datasets have begun to place useful 
constraints on global plate motions. It is clear, however, that efficient algorithms to 

handle the rapidly accumulating space-geodetic data, which are capable of treating 
the problems associated with each technique in a self-consistent manner, will be 

crucial for the detection and monitoring of temporal variations in plate motions. 
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TABLE 1. SLfc ANNUAL SOLUTION STATION SUMMARY 
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Epoch-Subset Results 


Epochs 

#Data 

Y 

°Y 

A 

s 2 

A 

1,2 

3 

-3.470 

0.450 

0.009 

22.540 

2,3 

11 

0.439 

1.090 

0.392 

3.028 

3,4 

8 

1.308 

0.673 

0.348 

1.303 

4,5 

11 

1.248 

0.372 

0.179 

0.770 

5,6 

9 

0.514 

0.501 

0.345 . 

0.726 

1,2,3 

15 

-0.226 

0.816 

0.368 

1.810 

2,3,4 

22 

0.905 

0.423 

0.315 

0.568 

3,4,5 

21 

1.275 

0.270 

0.279 

0.260 

4,5,6 

21 

0.969 

0.258 

0.277 

0.241 

1,2, 3, 4 

27 

0.587 

0.389 

0.363 

0.418 

2, 3, 4, 5 

35 

1.184 

0.232 

0.287 

0.187 

3, 4, 5, 6 

32 

1.120 

0.194 

0.341 

0.110 

1,2,3, 4, 5 

41 

0.963 

0.218 

0.327 

0.145 

2,3, 4, 5,6 

46 

1.061 

0.171 

0.335 

0.087 

1,2, 3, 4, 5, 6 

52 

0.947 

0.160 

0.355 

0.072 


Table 3. Results from estimation of- rate-scaling parameter based on epoch-subsets of 
the data. Epoch 1 = 197S, epoch 6 = 1983, etc; £Data = number of relative motion vectors 

observed during interval; the other quantities are defined in the text; all stations used 
for each estimation. 



Table 4. Comparison of Talwani et al. (1971) 
anomaly timescales, with corresponding 


and Ness et al. (1980) marine magnetic 
estimate of rate-scaling parameter. 


Anomalv 

Talwaniim.vO 

Nessfm.v.) 

y(T/N) 

1 

0.00 

0.00 

1.000 


0.69 

0.72 

0.958 

1' 

0.89 

0.91 

0.978 


0.95 

0.97 

0.979 

2 

1.71 

1.66 

1.030 


1.86 

1.87 

0.995 

2‘ 

2.43 

2.47 

0.934 


2.84 

2.91 

0.976 


2.94 

2.95 

0.997 


3.04 

3.07 

0.992 


3.10 

3.17 

0.978 


3.32 

3.40 

0.976 

3 

3.78 

3.86 

0.979 


3.88 

3.98 

0.975 


4.01 

4.12 

0.973 


4.17 

4.26 

0.979 


4.31 

4.41 

0.977 


4.41 

4.49 

0.982 


4.48 

4.59 

0.976 


4.66 

4.79 

am 
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Figure Captions 

Figure 1. Nominal SLR station positions used to estimate rate-scaling parameter. 
Station locations corresponding to site number are listed in Table 1. 

Figure 2. Data, consisting of apparent tangential vector motions of network sites 
positions between successive epochs with associated tangential one-sigma 
(unsealed) error ellipse, used to estimate rate-scaling parameter. Also listed in 

Table 2. 

Figure 3. Estimate of rate-scaling parameter as a function of total number of motion 

vectors utilized in epoch-subset trials. Dashed line is the RM2 reference estimate. 

Error bars are one-sigma (scaled by sampling variance). Values listed in Table 3. 

Figure 4. Estimate of rate-scaling parameter as a function of the temporal averaging 
interval. RM2 is the reference estimate with 1 m.y. error bars. One- and 

five-year SLR-derived estimates have errors as in Figure 3. Values listed in Table 
3. NLC-80 are estimates derived from the comparison of the Talwani et al. (197 1) 
marine magnetic anomaly timescale, upon which the RM2 rates are based, and 

Ness et al. (1980) timescale, as listed in Table 4. Variability of NLC-80 about RM2 is 
well within the five-year SLR-derived one-sigma level of confidence. 
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